Hydrodynamics of probabilistic ballistic annihilation 



Francois Coppex, 1 Michel Droz, 2 and Emmanuel Trizac 3 

department of Theoretical Physics, University of Geneve, CH-1211 Geneve 4, Switzerland 
2 Department of Physics, University of Geneve, CH-1211 Geneve 4, Switzerland 
3 Laboratoire de Physique Theorique (UMR 8627 du CNRS), 
Batiment 210, Universite de Paris-Sud, 91405 Orsay, France 

We consider a dilute gas of hard spheres in dimension d > 2 that upon collision either annihilate 
with probability p or undergo an elastic scattering with probability 1 — p. For such a system neither 
mass, momentum, nor kinetic energy are conserved quantities. We establish the hydrodynamic 
equations from the Boltzmann equation description. Within the Chapman-Enskog scheme, we 
determine the transport coefficients up to Navier-Stokes order, and give the closed set of equations 
for the hydrodynamic fields chosen for the above coarse grained description (density, momentum 
and kinetic temperature). Linear stability analysis is performed, and the conditions of stability for 
the local fields are discussed. 

PACS numbers: 05.20.Dd,47.20.-k,45.05.+x,82.20.Nk 

I. INTRODUCTION 

The hydrodynamic description of a low density gas of elastic hard spheres supported by an underlying kinetic 
theory attracted a lot of attention already more than 40 years ago 0, 0, @ ■ It nas now become a well established 
description ,4] . A key ingredient in the hydrodynamic approach is the existence of collisional invariants (quantities 
conserved by the -instantaneous- collisions). The question of the relevance of a coarse-grained hydrodynamic approach 
is therefore more problematic when the kinetic energy is no longer a collisional invariant [j| . This is the case of rapid 
granular flows (that may be modeled by inelastic hard spheres in a first approach), where the hydrodynamic picture, 
despite including a hydrodynamic field associated with the kinetic energy density, is nevertheless reliable (see, e.g., 
@,Q>H>IS an d Dufty for a review |10|). It seems natural to test hydrodynamic- like approaches further and in more 
extreme conditions, and investigate a system where particles react so that there exist no collisional invariants. The 
present article, focusing on the derivation of the hydrodynamic description for such a system, is a first step in this 
direction. 

The system we consider is made of an assembly of hard spheres that move ballistically between collisions. Whenever 
two particles meet, they either annihilate with probability p, or undergo an elastic collision with probability (1 — p). 
The model of probabilistic ballistic annihilation in one dimension for bimodal discrete initial velocity distributions 
was introduced in whereas for higher dimensions and arbitrary continuous initial velocity distributions it was 
considered in [l^. When p = 1, we recover the annihilation model originally defined by Elskens and Frisch [l3j . 
that has attracted some attention since 0, 0, 0, 0, 0, 0, |2(], l2lll In one dimension (again for p = I), the 
problem is well understood for discrete initial velocity distributions |l5l ITt| . On the contrary, higher dimensions 
introduce complications that make the problem much more difficult to treat p^l l2(ij |. Only a few specific initial 
velocity conditions lead to systems that are tractable using the standard tools of kinetic theory |2l| . 

Our starting point will be the Boltzmann equation, which describes correctly the low density limit of granular gases 
(see j^l and [23, 01 f° r the elastic case) . For annihilation dynamics, the ratio of particle diameter to mean free path 
vanishes in the long-time limit, such that for d > I the Boltzmann equation is valid at least at late times [2(],Ej. For 
such a dynamics, none of the standard hydrodynamic fields (i.e., mass, momentum, and energy) is associated with a 
conserved quantity. There are therefore three nonzero decay rates, one for each field. Numerical evidences show that 
in the long-time limit, a non-Maxwellian scaling solution for the homogeneous system appears (homogeneous cooling 
state (HCS) which also exists for inelastic hard spheres 25]). Nothing is known about the stability of the latter 

solution, the only existing result being that in one dimension, with a bimodal initial velocity distribution, clusters 
of particles are spontaneously formed by the dynamics |16|. In view of this situation we develop a hydrodynamic 
description for probabilistic ballistic annihilation. The limiting case of vanishing annihilation probability p — > gives 
the known results for elastic dilute gases 26], whereas the other limit p — *■ 1 yields the case of pure annihilation. 

In order to derive the hydrodynamic equations, we make use of the Chapman-Enskog method. We thus consider (at 
least) two distinct timescales in the system. The microscopic timescale is characterized by the average collision time 
and the corresponding length scale defined by the mean free path. The macroscopic timescale is characterized by the 
typical time of evolution of the hydrodynamic fields and their inhomogeneities. The fact that those two timescales are 
very different implies that on the microscopic timescale the hydrodynamic fields vary only slightly. Therefore they are 
on such length- and timescales only very weakly inhomogeneous. Combined with the existence of a normal solution 
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for the velocity distribution function (i.e., a solution such that all time dependence may be expressed through the 
hydrodynamic fields), this allows for a series expansion in orders of the gradients, i.e., application of the Chapman- 
Enskog method. The knowledge of the hydrodynamic equations thus obtained to first order allows us to perform a 
stability analysis. Taking the HCS as a reference state, we study the corresponding small spatial deviations of the 
hydrodynamic fields. 

The paper is organized as follows. In Sec. [H] we present the Boltzmann equation for probabilistic ballistic annihi- 
lation, and establish the subsequent balance equations. Sec. IIIII is devoted to the Chapman-Enskog solution of the 
balance equations. For this purpose we consider an expansion of the latter equations in a small formal parameter. 
The solution to zeroth order provides the hydrodynamic fields of the HCS. Assuming small spatial inhomogeneities, 
we make use of an explicit normal solution for the velocity distribution function to first order. This allows us to 
establish the expression for the transport coefficients and for the decay rates to first order, and thus the closed set 
of equations for the hydrodynamic fields. The technical aspects of the derivations are presented in the Appendixes 
while our main results are gathered in Eqs. I|47[) . In Sec. IIVI we study the linear stability of those equations around 
the HCS. Finally, Sec. Ivl contains the discussion of the results and our conclusions. Since from the point of view 
of dissipation probabilistic ballistic annihilation shares some features with granular gases, making several parallels 
between those two systems will prove to be instructive. 



II. BOLTZMANN AND BALANCE EQUATIONS 

The Boltzmann equation for the one particle distribution function /(r, v; t) of a low-density system of hard spheres 
annihilating with probability p is given by 

(ft + vx ■ V)/(r, vi; t) = pJalf, /] + (!- p)J c [f, /], (1) 
where the annihilation operator J is defined by |20j| 

J a [f,g] = -a d - 1 /3 1 [ dv a « ia /(r,v 2 ;t)0(r,vi;t) (2) 

and the elastic collision operator J c is defined by (2^, H3, H3 

Jc[f, g] = o d ~ x f dv 2 fda(a- Vi 2 )6(-o- ■ v 12 )(fc- 1 - l) 5 (r, v i; t)f(v, v 2 ; t). (3) 

In the above expressions, a is the diameter of the particles, vi 2 = vi — v 2 the relative velocity, vyi — |vi 2 |, 9 is the 
Heaviside distribution, <x a unit vector joining the centers of two particles at collision and the corresponding integral 
is running over the solid an gle, (3\ = 7r^ _1 ^ 2 /r[(<i+ 1)/2] where T is the gamma function, and b^ 1 an operator acting 
on the velocities as follows |28|: 

6 _1 vi2 = via - 2(vi2 • (4) 



b 1 v 1 = v x - (v 12 ■ er)S\ (5) 

Since J c describes elastic collisions, this operator conserves the mass, momentum, and energy. On the other hand, J a 
describes the annihilation process and thus none of the previous quantities are conserved. 

In order to write hydrodynamic equations, we need to define the following local hydrodynamic fields: 

n(r,t) = / cZv/(r,v;t), (6a) 



u(r,i) = — [ dvv/(r,v;t), (6b) 
n(v, t) J Md 

= < Tur A I dvV2 /(r,v;t). (6c) 
n(v,t)k B d J R d 



where n(r,i), u(r, t), and T(r,t) are the local number density, velocity, and temperature, respectively. The definition 
of the temperature follows from the principle of equipartition of energy. In Eq. (|6c|l . ks is the Boltzmann constant 
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and V = v u(r. t) is the deviation from the mean flow velocity. The balance equations follow from integrating the 
moments 1, mv, and mv 2 /2 with weight given by the Boltzmann equation Q). We thus obtain 

d t n + Vi(rati) = -pu [/,/], (7a) 

d t u i + —V j P ij +u j V j u i = -p-u[f,V i f), i = l,...,d, (7b) 
mn n 

2 T m 

where we have summation over repeated indices, u = [u\, . . . , Ud), and 

u[f,g] = <J d ~ 1 lh f dvidv 2 \vi2\g(vi,vi;t)f(r 2 ,v 2 ;t). (8) 

JM 2d 

In the balance equations (0, the pressure tensor P,j and heat-flux qi are defined by 

P ij (r,t)=m[ dvV i V j f(r,v;t)= [ dw f(r, v; t)A,-(V) + (9) 

<?<(r,f) = / dv^(V)/(r,v;i), (10) 

•/R d 



where /3 = l/{k B T) and 



£> ij (V) = m(F i F i -^-^ i ) ! (ill 



^(V)=(yy 2 -^ B Tj^. (12) 

One sees from Eqs. (|7J that when the annihilation probability p — > 0, all quantities are conserved. In addition, the 
long time solution of the system in this limit is given by the Maxwell distribution [l^ . 

III. CHAPMAN-ENSKOG SOLUTION 

In order to solve Eqs. Q, it is necessary to obtain a closed set of equations for the hydrodynamic fields. This can 
be done using the Chapman-Enskog method, by expressing the functional dependence of the pressure tensor Rj and 
of the heat flux qi in terms of the hydrodynamic fields. Note that other routes have been dev elop ed as well [271 129| . 
A thorough comparison of the different approaches seems, however, to be a difficult attempt 29] . In order to apply 
the Chapman-Enskog method, it is necessary to make two assumptions. The first one is that all temporal and spatial 
dependence of the distribution function /(r, v;t) may be expressed as a functional dependence on the hydrodynamic 
fields: 

/(r, v; t) = / [v, n(r, t), u(r, t), T(r, t)] . (13) 

What is the physical justification for the existence of such a normal solution? Suppose that the variations of the 
hydrodynamic fields are small on the scale of the mean free path, for example (na d ^ 1 )^ 1 \ V lnn| <§; 1. Therefore, 
to first order the functional dependence of the distribution function may be made local in the hydrodynamic fields, 
leading to the normal solution written above. Note that none of the hydrodynamic fields is associated with a conserved 
quantity. The theoretical question that arises is to know if the new timescales thereby introduced by the cooling rates 
are shorter than what is allowed for the existence of a normal solution ^3 ■ F° r sufficiently small p this should be 
the case. However, in the related context of granular gases, this point is not yet quantitatively clarified and is still 
subject to discussions 0, |3(], UD- The justification of the normal solution m ay b e done a posteriori by studying 
the relevance of the results through the appearance of the HCS, for example ^3>|2^|. The second assumption is 
based on the existence of (at least) two distinct timescales. The microscopic timescale is characterized by the average 
collision time and the spatial length which is defined by the corresponding mean free path. On the other hand, the 
macroscopic timescale is defined by a typical timescale describing the evolution of the hydrodynamic fields and their 
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inhomogeneities. The difference in those two timescaies implies that on the microscopic timescale the hydrodynamic 
fields vary only very slightly. Thus, those fields are on such time and space scales only very weakly inhomogeneous. 
This allows for a series expansion in orders of the gradients of the fields: 



/ = / (0) + A/ 



(i) 



+ A 2 / 



(2) 



(14) 



where each power of the formal small parameter A means a given order in a spatial gradient. The formal parameter 
A may be seen as the ratio of the mean free path to the wavelength of the variation of the hydrodynamic fields. This 
shows again the idea of the separation of both microscopic and macroscopic time and length scales. The Chapman- 
Enskog method assumes the existence of a timescale hierarchy, and thus of a time derivative hierarchy as well: 



at 



Q(.0) q(1) 



+ A : 



at 



(15) 



where a given order in the temporal hierarchy (|15|) corresponds to the same order in the spatial hierarchy (|14fl . One 
thus concludes that the higher the order of the spatial gradient, the slower the corresponding temporal variation. 
Inserting expansions H14[) and (|15fl in the Boltzmann equation one obtains 



E Afe V +Vl - v IE^ 



(0 



,k>Q 



■ pJa 



l>0 



+ (l-p)Jc 



(0 



z>o 



(16) 



Collecting the terms of a given order in A and solving the equations order by order allows us to build the Chapman- 
Enskog solution. 



A. Zeroth order 



To zeroth order in the gradients, Eq. I|16|l gives 

5 f (0) / (0) =pJa[/ (0) ,/ (0) ] + (l-p)J C [/ (0) ,/ (0) ]. (17) 

This equation has a solution, describing the HCS, and which obeys the scaling relation 

/(°)(r,v;t) = 7(c). (18) 

The approximate expression for /(c) was established in ^| and is recalled by Eq. 14011 . In Eq. I|18|) . vt — [2/(/3m)] 1 / 2 
is the time dependent thermal velocity, and c = V/vt, V = v — u. The existence of a scaling solution of the 
form H18JI seems to be a general feature that is confirmed numerically (direct Monte-Carlo simulations or molecular 
dynamics) not only for ballistic annihilation j20J or granular gases |32j , but for the dynamics of ballistic aggregation 
as well [3l EH- 

The function /(°) is isotropic. Thus to this order the pressure tensor © becomes = p^Sij, where — nksT 
is the hydrostatic pressure, and the heat flux 11UI) becomes q^ - 1 = 0. 

(19a) 
(19b) 
(19c) 



d t n = 


~pn^\ 




d t Ui = 




i = l,...,d, 


d f T = 


-PT^\ 





where the decay rates are 



<r(0) . 

Sn 


^[/ (0) ,/ (0) ], 
n 




t(0) _ 


— ^[/ (0) ,^/ (0) ], 
uvt 


i = l,. 


AO) _ 
?T — 


, m T ,-[/ (0) ^ 2 / (0) ] 


n 



(20a) 

,d, (20b) 



n (20c) 

It 

For antisymmetry reasons, one sees from Eq. (I20b|l that = 0. The two other decay rates are given later on by 
Eqs. gSl- 
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B. First order 



To first order in the gradients, the Boltzmann equation (|16fl reads 

[9 f (0) + J]/ (1) = -[a t (1) + v 1 .v]/(°), 

where 

Jf^=pL a [f^,f^] + {l-p)L c {f^,f^], 

with 

La[f {0 \f {1) } = -Ja[f (0 \f {1) }~Ja[f (1 \f { % 



(21) 
(22) 
(23) 



i C [/ (0) ,/ (1) ]=-J C [/ (0) ,/ (1) ]-Jc[/ (1) ,/ (0) ]. 

The balance equations J7J) to first order become 

d[ l) Ul + — V,(nT) + UjWjUt = -put$, 
ran 



i = l,...,d, 



d[ 1] T + u t V,T + -TViUi = -pT^\ 



where the decay rates are given by 



4 1} - — ^[/ (0 U/ (1) ] + — <*\f {l W% i = i,...,d, 



uvt 
rn 



= -- w [/(°^«] + ^ / [/'°^¥ 1) ] + 1 ^/[/ (1 ^ 2 / (0) ' 

n tikbI a ukbi a 



(24) 

(25a) 
(25b) 

(25c) 

(26a) 
(26b) 
(26c) 



By definition we know that /W must be of first order in the gradients of the hydrodynamic fields; therefore for a 
low density gas 



/C 1 ) = AiVi In T + BiV, In n + CyVj-? 



(27) 



The coefficients Ai, Bi, and C y - depend on the fields n, V, and T. Inserting Eq. I|27|l in Eq. I|21|) and making use of 
Eqs. (|13|1 . (|18|l . and (|19|l . one obtains the following set of equations for Ai, Bi, and C y - (see Appendix 0): 



{-p [^ 0) T9 T + ^nd n + ±4° } ] + ( J — P n)} Ai - p^Bi = A h 
[-p [^ ] Td T + ^ } nd n + ^ 0) ] + (J — pV)}Bi - p^Ai = Bi, 

{ -p [4 0) ra T + £i 0) na„] + ( j - pn) } Cy = c y , 



where 



a- Vi 9 iv m kBTdfl0) 

Vif{0) _kBTdf<» 
m dVi ' 
d „ r 1 d 

ddV k 



D 



_£_rv./(o)i _ i_£_rvLfC°)i5-- 



and fHs a linear operator defined by 



<9T 



(28a) 
(28b) 
(28c) 

(29a) 

(29b) 
(29c) 

(30) 



where g is either At, Bi, or Cy, and the functionals £« , £««"') an( i are obtained from Eqs. I|26|) upon replacing 
by g. It is possible to show that from Eqs (|29|l the solubility conditions ensuring the existence of the functions 
Ai, Bi, and are satisfied (see, e.g., H^|). 



(i) 
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C. Navier-Stokes transport coefficients 



The hydro-dynamic description of the flow requires the knowledge of transport coefficients. The concern of the 
present section is to determine the form and coefficients of the constitutive equations. This can thus be achieved by 
linking those macroscopic transport coefficients with their microscopic definition. Using a first order Sonine polynomial 
expansion, it is then possible to find explicitly the transport coefficients to first order. This will allow us to express 
the functions Ai, Bi, and Cy in terms of the transport coefficients, thus determining the distribution function Z*- 1 -*. 

The pressure tensor may be put in the form 



Pij(r, t) = p^S i:j - tj + VjUi - ~S i:i V* 



(31) 



where p^ = nk B T is the ideal gas pressure, and 77 is the shear viscosity. For a low density gas, the bulk viscosity 
£ vanishes; therefore the last term in the pressure tensor may be neglected [Tol l35l l36| . Fourier's linear law for heat 
conduction is 



-KViT — /iV^n, 



(32) 



where k is the thermal conductivity and \i a transport coefficient that has no analogue in the elastic case. A similar 
quantity appears for granular gases, which again is non- vanishing in the inelastic case only pri l37j. 
The identification of Eq. (|31|l with Eq. 10 using the result of the first order calculation yields 



Similarly, the identification of Eq. (|32(l with Eq. Ill 01) using the first order calculation leads to 



dvSi(V)/ 



(i) 



The main steps of the calculation are shown in Appendix [51 and the result is 



= JL 

Vo 

K 
Kq 

Tk 



1 



2 Kn 



V 



-(2a 2 



2v*. 



3*4 0> 



2 P ^ 



PC 



(Oh 



* 1 

K + 



d-1 



-a 2 



where ai is the kurtosis of the distribution 

a-2 = 

and 



d(d + 2) v: 



1- / dV/<°>(V)-l, 



d(d + 2) k B 
2(d — 1) m 



(33) 

(34) 

(35a) 
(35b) 
(35c) 

(36) 
(37) 



Vo 



d + 2 T(d/2) ^/mk B T 

~^~ n {d-l)/2 fjd-l ' 



(38) 



are the thermal conductivity and shear viscosity coefficients for hard spheres, respectively j3^|. £n°^* = £n° /vq and 
4°)*=4 0) MarethedimeE 
v* 1/*, and v* are given by 



£p * = $ / v o are the dimensionless decay rates, where isq = p^ /t)q, with p^ — nk B T . The dimensionless coefficients 



v* = 



1 f Kd dVSi(V)JAi 1 J^dVSiWSlAi 



1 hdVSiWW 



is SgtdVSiWAi ' 

1 L d dySi(v)m 



iso fjpdVSiMBi F iso $ Rd dVS % {V)B % ' 
1 L d dV D tJ (V)JCy 1 / dV D iS (VjnCi. 



iso / B< ,dV£>tf(V)Cy *V J^dVD^iV)^ 



(39a) 
(39b) 
(39c) 
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It must be emphasized that the above results are still exact within the Chapman-Enskog expansion framework. 
However, the relations (j39|l and the decay rates (|20[1 cannot be evaluated analytically without approximations. For 
this purpose, we first consider the Sonine expansion for /(°). It was shown that to first non-Gaussian contribution in 
Sonine polynomials the distribution /(°) reads [lj] 



/(°>(V) = 4a4 



where 



— 1 



«2 



IV 4 d + 2V 2 d(d + 2) 



2 v\ 



(40) 



M 



T d/2 



-W4 



(41) 



is the Maxwellian and 



«2 



Z-2V2 



Ad + 6 - y/2 + ^8V2(d - 1) 



(42) 



The coefficient a 2 was shown to be in very good agreement with direct Monte Carlo simulations [T^ ■ The relation l|4U|) 
allows us to compute the decay rates (see Appendix 0: 



(0)* _ 



JO)* 
ST 



d + 2 
4 

d+2 



1 + a 2 



8d+ll 
16 - 



Next, we retain only the first order in a Sonine polynomial expansion applied to *4., 0, and C. We thus have 

A(V) = ai M(V)S(V), 
B(V) = 6iM(V)S(V), 
C(V) = c .M(V)D(V), 



(43a) 
(43b) 



(44a) 
(44b) 
(44c) 



where a%, bi, and cq are the coefficients of the development. This allows us to compute the relations l|39|l . For this 
purpose, as already shown the probabilistic collision operator J given by Eq. (1221) can be split into the sum of an 
annihilation operator and of a collision operator. Each contribution may thus be treated separately. Therefore we 
make use of previous calculations for the collision process |39| . The calculations for the annihilation operator are 
shown in Appendix [D] and the final results read 



1 



V K = V U = P 



32d 
+(l-p) 



16 + 27d + 8d 2 
d-1 



a 2 - 



2880 + 1544d - 2658d 2 - 1539d 3 - 200d 4 



32d(d+2) 



1 + a 2 — 
32 



1 

^ = P 8d 



3 + 6d + 2d 2 - a 2 



278 + 375d + 96d 2 + 2d 3 



32(d + 2) 



+ (l-p)(l-a 2 - 



(45a) 
(45b) 



One may check that these expressions approach the known elastic values when p — > 0] 39] . The transport coefficients 
are thus found from Eqs. (JHSJ using Eqs. J37J|, JSHJ, (gSJl, and |g5J|. 

In order to establish the decay rates to first order, one needs the distribution /W (see Appendix 0: 



/W(r,V;t) = -^M(V) 
n 



2m 
d + 2 



(46) 



D. Hydrodynamic equations 



The pressure tensor and the heat flux defined by Eqs. H31JI and H32JI ■ respectively, are of order 1 in the gradients. Thus 
their insertion in the balance equations yields contributions of order 2 in the gradients. Consequently there are 
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second order terms (so called Burnett order) that contribute to the first order (so called Navier-Stokes order) transport 
coefficients, and the knowledge of the distribution is thus necessary. Indeed, use was made of the zeroth order 
relations — p^°'Sij and qi — to establish the balance equation for energy (|25c|l . However, it was shown in the 
framework of the weakly inelastic gas - and consequently for an elastic gas - that those Burnett contributions were 
three orders of magnitude smaller than the Navier-Stokes contributions 26] . For the sake of simplicity, we will here 
neglect those second order contributions. For small annihilation probabilities p, this approximation is thus likely to 
be justified. However, we have a priori no control on the error made when the annihilation probability p is close to 
unity. 

The hydrodynamic Navier-Stokes equations are given by 



d t n + V i (nu i ) = -pn[^ + &% 



1 



dtUi -\ VjPy + uNjUi = -pv T £,u- 

mn 

d t T + u % V,T 



i = 1, 



nksd 



(47a) 
(47b) 

(47c) 



where the decay rates £n and £y/ are given by Eqs. I|43a|) and l|43b|l . respectively. and qj are given by Eqs. 

with C = 0, and !i'S3 respectively. The rates £,n \ , and $p may be calculated using their definition l|25[l and the 
distribution l|46|). We find (see Appendix IF|l : 



= 



AD 
St 



T 



0, 



(48a) 
(48b) 

(48c) 



where 



(d + 2) 2 
32(d- 1) 



1 + a 2 - 



-86- 101d+32d 2 



32(d + 2) 



28d 



'41 



(49) 



We thus have a closed set of equations for the hydrodynamic fields to the Navier-Stokes order. 



IV. STABILITY ANALYSIS 



The hydrodynamic Eqs. I|47|) form a set of first order nonlinear partial differential equations that cannot be solved 
analytically in general. However, their linear stability analysis allows us to answer the question of formation of 
inhomogeneities. The scope of the present study is to find under which conditions the homogeneous solution to zeroth 
order, i.e., the HCS, is unstable under spatial perturbations. To this end we consider a small deviation from the HCS 
and the linearization of Eqs. (|47|l in the latter perturbation. Equations l|19|) give the time evolution of the HCS, which 
is found to be 

(50a) 

It 

(50b) 

where the decay exponents are j n = ^n\o)to, 7t = £j?^(0)ioj and the relaxation time i^ 1 = £n\o) + £^(0)/2. The 
subscript H denotes a quantity evaluated in the homogeneous state. The density and temperature fields of the HCS 
are thus decreasing monotonically in time, with exponents that depend on the annihilation probability through the 
kurtosis of the velocity distribution. The explicit expression of the decay exponents may be obtained straightforwardly 
using Eqs. |@3J|. 

The linearization procedure used here follows the same route as the method used for granular gases [2(J . We define 
the deviations of the hydrodynamic fields from the HCS by 



n H (t) = n [l +P— 



T H (t) = T [l+p- 



Sy(r,t) = y{r,t) - y H (t), 



(51) 
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where y = {n, u, T}. Inserting the form (|51|l in Eqs. I|48(l yields differential equations with time-dependent coefficients. 
In order to obtain coefficients that do not depend on time, it is necessary to introduce the new dimensionless space 
and time scales defined by 



( j fji 

1 = ? aH{t) ikMT) v > 



as well as the dimensionless Fourier fields 



where 



Pk(r) 
w k (r) 

<5yk(-r) 



dsvoH(s), 
u h (t) ' 

/ TO 
V k B T H (T) 

5T k (r) 
Th{t) ' 



5uk(r), 



die 



-ik-1 



52/(1, r) 



(52a) 
(52b) 

(53a) 
(53b) 
(53c) 

(54) 



From Eq. (|52a[l , it appears that lengths are made dimensionless making use of the time dependent mean free path as 
a reference scale. Making use of Eqs. (|53|) and H52fl in Eqs. I|47|). the linearized hydrodynamic equations read 



d_ 

dr 



Mil 



(0)* 



Pk(r) +Ki 0) ^k(r) + i fc Wk|1 (r) = 0, 



^7 



(0), 



*4 0) * 



d + 2 
2(d- 1) 



- * ; 2 
-?7 K 



w k|1 



ik 



(1-K>*)pk(r) + (I-K>*)^k(r) 



= 0, 



« p€t 

OT 



(0)* 



1 * 7 2 

2 ' 



w kx (r) = 0, 



-*h 2 



c(r) 



d + 2 
2(d-l 



* 7 2 
/j A; 



Pk(T") + -jikw^ (t) = 0, 



(55a) 
(55b) 
(55c) 
(55d) 



where wt. and are the longitudinal and transverse part of the velocity vector defined by w k|| = (w k • e k )e k 
and w ki = w k — w k|| , where e k is the unit vector along the direction given by k. Eq. I|55c|l for the shear mode is 
decoupled from the other equations and can be integrated directly so that 



w kj _(r) = w k _ L (0)exp[sj.(p,fc)r], 



where 



s±(p,k) = p$ ] * - - 



r]*k 2 



(56) 



(57) 



The transversal velocity field w ki lies in the (d — 1) dimensional vector space that is orthogonal to the vector space 
generated by k, and therefore the mode s± identifies (d — 1) degenerated perpendicular modes. The longitudinal 
velocity field w k|| lies in the vector space of dimension I generated by k. Hence there are three hydrodynamic fields 
to be determined, namely the density p k , temperature # k , and longitudinal velocity field w k|| = iu k|! e k . The linear 
system thus reads 




(58) 



with the hydrodynamic matrix 



(0)* 



M 



-2p£ 

-tk(i- P Cn* 



d+2 



—ik 

d-l „* L.2 



(0)* 



-ih(l-pQK*) 



lik 



(59) 
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FIG. 1: Real part of the eigenvalues in dimensionless units for probabilistic ballistic annihilation with p = 0.1 and d — 3. The 
dispersion relation obtained from Eq. 157H is represented by a dashed line (labeled s±) whereas the three remaining relations 
obtained upon solving Eq. 1591 are represented by continuous lines (labeled sit). 

The corresponding eigenmodes are given by <p n (k) = exp[s„(p, fc)r], n = 1, . . . , 3, where s n (p, k) are the eigenvalues 
of M. Each of the three fields above is a linear combination of the eigenmodes; thus only the biggest real part of the 
eigenvalue s n (p, k) has to be taken into account to discuss the limit of marginal stability of the parallel mode of the 
velocity field. Figure. Q shows the real part of the eigenvalues for p — 0.1 and d — 3 (obtained numerically). 

One may identify three regions from the dispersion relations. We first define k± (dimensionless) by the condition 
Re[s±{k±,p)} = 0, i.e., 



(60) 



and by max^ Re[s||(fc||,p)] = (the expression for k\\ is too cumbersome to be given here); we have < k±. 
Figure [5] shows the dependence of k± and fen as a function of the annihilation probability p. Then for all k > k± all 
eigenvalues are negative and therefore, according to Eq. Ij56(l . correspond to linearly stable modes. For k G [feikfcjj 
only the shear mode Wi Ii of the velocity field is linearly unstable. In the case of granular gases in dimension 
larger than one this region exhibits velocity vortices (2^, l4pL El) , with a possible subsequent non- linear coupling to 

density inhomogeneities. From = CtV^o an( i Eq- l|52b|) one may integrate Eq. (|19c|) in order to find Th{t) = 
Tff(0) exp[-2pd 0) *T]. Then equating Eqs. 
and of Eq. I|53bj) for r = 0, one finds 



53b|l and (|56|l . making use of the latter expression for Th(t), of Eq. JSTJ), 



5u k± (T) = u ki (0)exp ( --r]*k 2 r 



(61) 



The exponential decay in the reduced variable r translates into a power-law-like decay in the original variable t [since 
the exponent k = k(t) depends itself on time]. Indeed, the integration of Eq. (|19c|) yields t = — ln[T#(t)/TH-(0)]/2£, 



(0)* 



which we replace in Eq. I)61|l and make use of the homogeneous solution Ty(t) given by Eq. 
obtain 



Q) in order to finally 



<>u u j/) = u k _(.t)l j 



(62) 
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FIG. 2: Wavenumber k± and fen in dimensionless units as a function of the annihilation probability p for d = 3. 

where ig = io/^-f/W is the dimensionless relaxation time. In the linear approximation the perturbation of the 
transversal velocity field therefore decays even if s±(k,p) > 0. The rescaled modes with k < kn are linearly unstable. 

However, a crucial point is that for any real (finite) system, the wave numbers are larger than 2n/L (assuming 
a cubic box of size L), which corresponds to a time dependent dimensionless wave number fc m j n = 27r/(Lner d ~ 1 ), 
which increases with time as 1/n. This lower cutoff therefore inevitably enters into the stable region k m [ n > kj_ , 
so that an instability may only be a transient effect. In other words, an unstable mode associated with a given 
value of k corresponds to a perturbation at a wavelength which increases with time in real space, and ultimately 
becomes larger than system size. However, at late times, the Knudsen number defined as the ratio of mean free 
path (which is proportional to k m i n ) over system size, becomes large, which should invalidate a Navier-Stokes-like 
description. Similarly, the present coarse-grained approach is a priori restricted to low enough values of k. Given that 
k± increases quite rapidly with p (see Fig.[2J|, the stable region k > k± might correspond to a "non-hydrodynamic" 
regime when p is larger than some (difficult to quantify) threshold. Conclusions concerning the stability of the system 
for such parameters rely on the validity of the hydrodynamic description (which could be tested by Monte Carlo or 
molecular dynamics simulations) which is beyond the scope of the present article. 

At this point, we conclude that the system may exhibit transient instabilities, but safe statements may only be 
made for very low values of p for which k±_ is low enough to guarantee that the hydrodynamic analysis holds. The 
stable region is then ultimately met irrespective of system size. 

With the above possible restrictions in mind, it is instructive to consider the counterpart of Fig. for "large" 
values of p (see Fig. For p > 0.893 . . ., we obtain the unphysical result that some eigenvalues increase and diverge 
upon increasing k. An a priori similar deficiency was reported for the inelastic Maxwell model where some transport 
coefficients — which are obtained exactly within the Chapman-Enskog method — diverge for strong dissipation 32]. 

V. CONCLUSION 

In this paper we construct a hydrodynamic description for probabilistic ballistic annihilation in arbitrary dimension 
d > 2, where none of the hydrodynamic fields can be associated with a conserved quantity. The motivation is not 
only to discuss the possibility of large scale instabilities in such a system, but also to provide the starting point 
for further (numerical) studies centered on the applicability of hydrodynamics to systems in which there are no 
collisional invariants. To this aim, we consider the low density and long-time regimes in order to make use of the 
Boltzmann equation with the homogeneous cooling state (HCS) as a reference state. The Chapman-Enskog method 
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FIG. 3: Real part of the eigenvalues in dimensionless units for probabilistic ballistic annihilation with p — 0.95 and d — 3. The 
figure caption is the same as for Fig. 

then allows us to build a systematic expansion in the gradients of the fields, with an associated timescale hierarchy. 
We consider only the first (Navier-Stokes) order in the gradients to build the hydrodynamic equations describing 
the dynamics of probabilistic ballistic annihilation. The transport coefficients and decay rates are established from 
the microscopic approach neglecting Burnett contributions and restricting ourselves to the first non-Gaussian term 
in a Soninc expansion. We then linearize the hydrodynamic equations around the HCS. The subsequent dispersion 
relations inform on the range of the perturbation's wavelength and time scales for which the system may exhibit 
density inhomogeneities. 

Interestingly the behavior of the dispersion relations and of the wave numbers k± and fc|| is qualitatively similar 
to its counterpart obtained for (inelastic) granular gases j2|| 0] . This leads us to conclude that some features of 
those models do not depend on the details of the dynamics, but rather on the parameter controlling the dissipation 
(referring to the existence of non conserved quantities) in the system, namely, p [or (1 — a 2 ) in the case of granular 
gases, where a is the restitution coefficient]. However, a specific feature of our model is that the mean free path 
increases rapidly with time. Consequently, even if the stability analysis leads us to the conclusion that this feature 
drives the system in a region where the homogeneous solution with a vanishing flow field is stable, the associated 
Knudsen numbers may be too large to validate our coarse-grained approach. At very small values of p, however, the 
stable region (k > k± with the notations of Sec. HV|) should be relevant, but then, the effects of transient instabilities 
in the case where the system is large enough to allow for fc m i n < k±_ or fc m j n < ku seem difficult to assess. 

Another point to emphasize is the amplitude of the dissipation in the system, which appears through the decay 
rates of the hydrodynamic fields. Again, there must be a clear separation between the macroscopic time scales 
described by those decay rates, and between the microscopic time scales. This separation of scales is required in 
the hydrodynamic approach in order to make use of the hydrodynamic fields n, u, and T that are associated with 
nonconserved quantities. The decay rates having the dimension of the inverse of a time, their inverse thus defines 
a time scale. If those decay rates increase, the associated time scales decay. In our case, we clearly introduce three 
such time scales that are supposed to be macroscopic; one for each nonconserved field. It is therefore required that 
the maximum of these decay rates defines a macroscopic time scale that is much bigger than the microscopic one. 
Nevertheless those decay rates increase as a function of the annihilation probability and hence the decrease of the 
associated time scale. One question that arises is to determine for which value of p the smallest time scale introduced 
by the decay rates is of the order of the microscopic time scale - which increases as a function of time because of the 
decreasing density of particles remaining in the system. When this is the case, the hydrodynamic description becomes 
irrelevant and one may not make use of the fields n, v, and T any more. As the parameter p controls the dissipation 
in the system, the question at hand here -left for future work- is reminiscent of the controversial issue of the validity 
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of hydrodynamics for granular gases with "low" coefficients of restitution. 
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APPENDIX A: INTEGRAL EQUATIONS FOR A, Bi, AND &, TO FIRST ORDER 

Using the normal solution (|13|) . the scaling form l|18[l. and the balance equations l|19|) one may rewrite the right 
hand side of Eq. (|21[1 in such a way that 



where Vi — — u,, 



and 



[<9 t (0) + J]/ (1) = AiVi InT + BiVi Inn + OyV^ +pSlfW, (Al) 



n/ w = /«»#) _ + ^r#>, (A2) 



4 = ^. ( A3a ) 

01 m oVi 

^ - + (A3C) 

The velocity dependence of /^°^ occurs only through V/vt- Because of the normalization the temperature dependence 
of the function is of the form T- d / 2 J {0) (U/T 1 / 2 ), and one obtains 

The insertion of Eq. (|A4() in Eqs. I|A3|I yields the relations lf2"§l) . 

Using the scaling form (|18|l and by definition of the decay rates <|2(J|) one has 

e<?> ~ 4 0) ~ "TV"- (as) 

This yields the relations 

t J §F = 5* (A6) 

and 

where = , Cy''}- Using the form (|27(l in the left hand side of the Boltzmann equation 1|21|) and making use 
of the relations l|A6(l and (|A7|) one obtains after some calculations the relations l|28|) . 
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APPENDIX B: EQUATIONS FOR THE TRANSPORT COEFFICIENTS 



As we will apply a Sonine expansion, the symmetry properties of «4.(V) and B(V) are the same as those of S(V), 
whereas the properties of C(V) are the same as those of D(V). Thus the insertion of Eq. I|27|) in Eq. (|33[1 yields 



= [ dvD«(V)C«CV)V fc ui. 
Jr* 

The identification of Eqs. (|B 1 1) and 131|1 yields (see e.g. |3fij ) 



(Bl) 



(B2) 



Integrating Eq. (|28c() over V in R d with weight — l/[(d — l)(d + 2)]Dy(V) and making use of Eq. (|B2|I one obtains 



p^Td T -p^nd n 



(d-i)(d+2)y B d 



dVD tJ (V)C K 



(B3) 



Functional dependence analysis shows that nd n r) — and Tdri] — r)/2. Using the definitions (|1 lfl for Dy(V) and (|29c(l 
for Cjj (V), it is possible to compute the right hand side of l|B3|) which gives 



1 



1 



dVmV 2 f(°l 



(B4) 



Using the hydrostatic pressure p(°) = nksT with the definition l |5c)l for the temperature, and dividing Eq. (|B4ll by rjo 
we finally obtain Eq. I|35al) . 

The insertion of Eq. JUJ in Eq. (PU gives 



q t = f dV Si(V)A k (V)V k lnT + f dV S,(V)6 t (V)V fc In, 

JR d JR d 

The identification of Eqs. I|B5|) and yields 

« J JR d 

H = ~ [ dV Si(V)Bi(V). 
dn J R d 



(B5) 

(B6a) 
(B6b) 



The fact that /i ^ is due to the annihilation process. Integrating Eqs. (|29a|) and 129b(l over V on K d with 
weight —Si(V)/d and making use of Eq. I|B6[) . then making use of T8t(Tk) — 3Tk/2, Tdrinfi) = 3n/i/2, and 
nd n {Tn) = nd n (nfj,) = obtained from functional dependence analysis, it follows 



1 



M = 



Using Eqs. I|29a|) . 129bl) . l(T2^l . and (|4*2|) one may calculate the integrals appearing in the right hand side of Eqs. (|B7() 



p£, ( t ] Tk- \ [ dVSi(V)Bi(y) 

« JM d 



1 



(B7a) 
(B7b) 



-L[ dv Si (v)Mv) = ~^(^2 i) 

dl JR d 



2 m/3 
d+2 1 



*()*() i ;-'/,! 

The insertion of Eqs. ifBSf in (fB7)l yields Eqs. (|35bjl and l|35c[) . 



2a 2 . 



(B8a) 
(B8b) 
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APPENDIX C: EVALUATION OF £i 0) * AND ^ 0) * 

The decay rates (|20a() and l|20c[) may be computed using the definition (JSJ and Eqs. (|40l) and (|41|l . We first change 
variables to Ci = Vi/vx, i = 1,2, then to c i2 = ci — c 2 and C = (ci + c 2 )/2 in order to decouple the integrals. Next, 
the integrals being isotropic with a symmetric weight, only even powers of the components of C and c 12 will give 
nonzero contributions. Thus the terms (C-Ci 2 ) 2 in the integrals become C 2 c 2 2 /d. Finally, the resulting integrals may 
be computed using the following relation |2J|: if we define 

M °nv = \l dc 12 dCe-^l 2 e- 2C \ n 12 C^ (CI) 

J -BP* 



M np = (4 2 C p ) 



\ [ dc 12 dCe- c ^ 2 e- 2cQ c? 2 c4l + aJc 4 + i-< 

7T a J R 2d [ { 16 



d + 2 
2d 



r 2 r 2 

O C 12 



(d + 2)C 2 - d -±^c 2 2 + d -^l 



then 



M _ n(n - p)/2 T[(d + n)/2}T[(d + p)/2} 

np r(d/2) 2 



> (C2) 



(C3) 



M np 1 r a 2 



M° p 16d 



Equations (|C3|) and (|C4ll may be easily verified using the relation 



, X|X| 6 - a d+n)/2 r(d/2) ' (C5) 



for a £ R + . We thus obtain the decay rates to zeroth order |@SJ|. 



APPENDIX D: EVALUATION OF u*, v*, AND v* 

Using the first order Sonine expansion l|44|) . Eqs. i|39|l reduces to 



1 J Rj( iVD j3 (V)J^A;] 



1 / R „ dV Dij(V)Sl[MDij] 



v / dV Aj ( V)M (V) Ai (V) ' / R2 dV Ai ( V)M (V)Aj (V) ! 



1 r Rd dV5 i (V)J[7W5 J ] 



i / Rd dvsi(v)n[^; 



i/ ; R2 dv Si(v)M(v)Si(v) % j R2 dv Si(v)M(v)Si(v) ■ 

The denominators of Eqs. (|D1|I are straightforward to compute using the formula l|C5(l . We thus find 

P 2 



(d + 2)(d- l)w 

2m/3 3 
d(d + 2)uvq 



[ dV Dij (V)J[M D i:j } - p [ dV Aj (V)O [M A 

/ dVSi(V)J[MSi]-p f dVS t (V)n[MS t 

Jm d JR d 



(Dla) 
(Dlb) 

(D2a) 
(D2b) 



The collision operator J defined by Eq. Ill'2[l is made of the sum of an annihilation operator L a with weight p and 
of an elastic collisional operator L c with weight (1 — p). Using previous calculations for L c [39| . we thus obtain the 
elastic gas contributions proportional to (1 — p) in the right hand side of Eqs. I|45|) . 

The following computations are technically simple, but lengthy. We shall thus only give the main steps. The 
annihilation contributions, written v* a , i/* a , and v* a , are given by 



2 



(d + 2)(d- l)nu 
2m/3 3 



dV Dij(V)L a [MDi 



d(d + 2)nvo J R d 



dV Si(V)L a [MSi] 



(D3a) 
(D3b) 



16 



where L a is given by Eqs. (HJ and ©, and 



K? = -TJTWT^S / dV D^VMMD^l (D4a) 

' (d + 2)(d- l)nv J R d 



< a = < 



2m0 



3 



d(d + 2)nvo 
Using the relation 



J dV Si(V)Q[MSi\. (D4b) 



/ dv 1 Y(v 1 )L a [MX]=a d - 1 /3 1 f d Vl dv 2 |v 12 |/( )( Vl )7W(v 2 )X(v 2 ) [Y(vi) + Y(v 2 )] , (D5) 

where X and V are arbitrary functions, changing variables to = v^/wt for i = 1,2, then changing variables to 
C12 = Ci c 2 [in the following we adopt the notation ci 2 = (ci 2l , . . . , ci 2d )] and C = (ci + c 2 )/2 in order to decouple 
the integrals, replacing under the integral sign for symmetry reasons the relations (C • Ci 2 ) 2 by C 2 c 2 2 /d, and using 

/ dCdc 12J F(C)G(ci 2 )(C-c 12 ) 4 = / dCdc 12 F{C)G(c 12 )(^C i c\ 2 -2dCtci 2 \ (D6) 

where i and j can be chosen arbitrarily in the set {1, . . . , d} and F, G are arbitrary isotropic integrable functions, one 
obtains 



< = MA \^Y^) HlM+ < a '' (D7a) 



d{d- 

i r(d/2) 



where 



H k (a 2 ,d)= a iJ [ dCe~ 2c2 C l [ dc^e^/ 2 ^ 1 

+ V 7ij / dCe-^Vc? / dci a e-^/ a ^M ai> (D8) 

with tty and 7-y that are functions of d and a 2 , and being the sets of allowed values for the pairs defining 
the moments in the integrals (|D8|I . Expressions for ay, 7«, f2*, and fii? are too cumbersome to be given here. The 
integrals in the first sum of the right hand side of Eq. (|D8I) may be computed using the formula (|C5|) ■ The integrals 
in the second sum may be computed using the particular case i = j = k = lol the formula 



Mg = / dxL\x\ n e- ax2 x iXj x k xi = 6<°> 



&ijki + g(l - hjki) (SijSki + 5 lk 8 3 i + SaSjk) 



(D9) 



where 



,(«)_ d/2 3 (d + n)(d + rc + 2)r[(d + n)/2] 1 

4 d(d + 2) r(d/2) a (rf+«+4)/2' { - uw > 

Using Eqs. (|C5() and l|F5(l we find = v* a = 0. The calculation can thus be performed and we obtain the first 
terms in the right hand side of Eqs. l|4"H|l. 

APPENDIX E: THE DISTRIBUTION 

Using the form (|27|l for /W and the first order Sonine expansion (|44|) one has 

/( 1 )(r ) V;i)=M(V)[a 1 5 i (V)ViT + 6 1 S , i (V)V i n + Co Ai(V)V^ i ]. (El) 
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The coefficients a\, b\, and cq may be expressed as functions of the transport coefficients, thus determining /W 
Eq. I|ij2l) in which we insert the Sonine expansion (|44|l yields 



V = ~ C ° (d-l)(d + 2) / dV ^( V )MV)A,(V) ' E- ! 



where we have made use of the definitions (|ll|l and l|41|) . 

Proceeding in a similar way with Eq. (|B6(I and (|44|) it follows 



-ax -J- / rfVS,(V)W(V)S,(V) 



(E3) 



(E4) 



d + 2nk B 

= -01—5 53, E5 

/i = f ^dV Si(V)M(V)Si(V) (E6) 

2 TOp J 

Replacing in Eq. IjElf) the coefficients cq, a\, and 61 obtained from Eqs. i|E3() . l|E5fl. and (|E7|I . one obtains the 
distribution ggj. 



APPENDIX F: EVALUATION OF ^ AND $ 



The zeroth order and first order distributions J*- ) and being known, it is possible to compute the first order 
decay rates ()26[). The procedure is similar to the calculation of the zero order decay rates of Appendix |0 We first 
change variables to Cj = Vi/vt, i = 1,2, then to C12 = ci — C2 and C = (ci + C2)/2 so that 



r (^±i) ir d 4 



Yd- 1 



i n 



(Fl) 



where u[f,g] is defined by Eq. ©, (A, B) = {(1, 1), (V 2 , 1), (1, V, 2 ), (V 2l ,l), (1, V, = (V$ 15 . . . , V id ), and 



7i=/ dc 12 |c 12 |e- c -/ 2 / dCe- 2c2 ^( VT c 2 ) J B( V TC 1 ) fc 2 -^) Cll [l + a 2 5 2 (c 2 )] , (F2) 



/a = 



/ dc 12 \c 12 \e~ c ^ 2 / dCc- 2c2 A(v T c 2 )B(v TCl ) ( c uCl] - ~<%c 2> ) [l + a 2 S 2 (c 2 2 )] , (F3) 



In the above integrals, ci and C2 are expressed as functions of C and C12. Then, in order to compute those integrals 
one needs the following additional relations: 



M^ = [ d X \x\ n e- ax2 x lXj = M^6 ij , 



(F4) 



where 



M («) _ d/2 d + nT[(d + n)/2] 1 

2d r(d/2) a («i+«+2)/2- [ - r ° ) 

Assuming summation over repeated indices, it is easy to show that 

M^M^ = M^M^6 ij , (F6) 
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(F7) 



Using the definition i|D9|l one can show that 




d + 2 
3 



b {a) M {a,) 5,, 



(F8) 




jklm 



d + 2 
3 



b {a) b {a,) 5,, 



(F9) 



M, 



m klmn ~ 3 



(F10) 



For symmetry reasons, many of the terms in the integrals l|F2(l and i|F3() vanish upon integration. Nevertheless, the 
expressions are very cumbersome and the use of symbolic computation programs is appreciable |43|. Equations. I|C5|I . 
(|D9|I . (|F4|) . and (|F6fl ~ (|F9ft thus allow us after a lengthy but technically simple calculation to find the decay rates to 
first order Q26|). 
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